home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / determ.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  112 lines
  1. ;$Id: determ.pro,v 1.10 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       DETERM
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the determinant of an N by N array.
  11. ;
  12. ; CATEGORY:
  13. ;       Linear Algebra.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;       Result = DETERM(A)
  17. ;
  18. ; INPUTS:
  19. ;       A:      An N by N array of type: float, or double.
  20. ;
  21. ; KEYWORD PARAMETERS:
  22. ;       CHECK:  If set to a non-zero value, A is checked for singularity.
  23. ;               The determinant of a singular array is returned as zero if
  24. ;               this keyword is set. Run-time errors may result if A is
  25. ;               singular and this keyword is not set.
  26. ;
  27. ;       DOUBLE: If set to a non-zero value, computations are done in
  28. ;               double precision arithmetic.
  29. ;
  30. ;       ZERO:   Use this keyword to set the value of floating-point
  31. ;               zero. A floating-point zero on the main diagonal of
  32. ;               a triangular matrix results in a zero determinant.
  33. ;               For single-precision inputs, the default value is 
  34. ;               1.0e-6. For double-precision inputs, the default value 
  35. ;               is 1.0e-12.
  36. ;
  37. ; EXAMPLE:
  38. ;       Define an array (a).
  39. ;         a = [[ 2.0,  1.0,  1.0], $
  40. ;              [ 4.0, -6.0,  0.0], $
  41. ;              [-2.0,  7.0,  2.0]]
  42. ;       Compute the determinant.
  43. ;         result = determ(a)
  44. ;       Note:
  45. ;            See CRAMER.PRO, in the same directory as this file, for
  46. ;            an application of the determinant function.
  47. ;
  48. ; PROCEDURE:
  49. ;       LU decomposition is used to represent the input array in
  50. ;       triangular form. The determinant is computed as the product
  51. ;       of diagonal elements of the triangular form. Row interchanges
  52. ;       are tracked during the LU decomposition to ensure the correct   
  53. ;       sign, + or - .
  54. ;
  55. ; REFERENCE:
  56. ;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
  57. ;       Erwin Kreyszig
  58. ;       ISBN 0-471-55380-8
  59. ;
  60. ; MODIFICATION HISTORY:
  61. ;       Written by:  GGS, RSI, February 1994
  62. ;       Modified:    GGS, RSI, November 1994
  63. ;                    Added CHECK keyword to check for singular arrays.
  64. ;                    Changed NR_LUDCMP to LUDC.
  65. ;       Modified:    GGS, RSI, April 1996
  66. ;                    Modified keyword checking and use of double precision.
  67. ;-
  68.  
  69. FUNCTION Determ, A, Check = Check, Double = Double, Zero = Zero
  70.  
  71.   ON_ERROR, 2  ;Return to caller if error occurs.
  72.  
  73.   Type = SIZE(A)
  74.   if Type[1] ne Type[2] then $
  75.     MESSAGE, "Input array must be square."
  76.  
  77.   if Type[3] ne 4 and Type[3] ne 5 then $
  78.     MESSAGE, "Input array must be float or double."
  79.  
  80.   ;If the DOUBLE keyword is not set then the internal precision and
  81.   ;result are identical to the type of input.
  82.   if N_ELEMENTS(Double) eq 0 then Double = (Type[Type[0]+1] eq 5)
  83.  
  84.   if N_ELEMENTS(Zero) eq 0 and Double eq 0 then $
  85.     Zero = 1.0e-6  ;Single-precision zero.
  86.   if N_ELEMENTS(Zero) eq 0 and Double ne 0 then $
  87.     Zero = 1.0d-12 ;Double-precision zero.
  88.  
  89.   if keyword_set(Check) then $ ;Return a determinant of zero?
  90.     if COND(A, Double = Double) eq -1 then $
  91.       if Double eq 0 then RETURN, 0.0 else RETURN, 0.0d
  92.  
  93.   ALUd = A ;Make a copy of the array for its LU decomposition.
  94.  
  95.   ;Compute LU decomposition.
  96.   LUDC, ALUd, Index, Double = Double, Interchanges = Sign
  97.  
  98.   ;Are there any zeros on the main diagonal?
  99.   ii = WHERE( ABS( ALUd[LINDGEN(Type[1])*(Type[1]+1)] ) le Zero, Cnt)
  100.  
  101.   Det = 1 ;Initialize determinant.
  102.  
  103.   if Cnt ne 0 then $ ;A zero on the main diagonal results in a zero determ.
  104.     if Double eq 0 then RETURN, 0.0 else RETURN, 0.0d $ 
  105.   else begin
  106.     FOR k = 0, Type[1]-1 do $
  107.       Det = Det * ALUd[k,k]
  108.     RETURN, Sign * Det
  109.   endelse
  110.  
  111. END
  112.